import Pkg

Pkg.activate("./")

using DistributedGISA, DelimitedFiles

function main(d::Int)
# parameters
Da1 = 1.0
Da2 = 2.0
dt = 1.0

# Reactor 1
function f1(x::AbstractArray, u::AbstractArray)
    x1 = x[1] + dt*(-x[1] - Da1*x[1] + u[1])
    x2 = x[2] + dt*(Da1*x[1] - x[2] - Da2*x[2]*x[2])
    return [x1 x2]
end

# Reactor 2
function f2(x::AbstractArray, u::AbstractArray)
    x1 = x[1] + dt*(-x[1] - Da1*x[1] + u[1])
    x2 = x[2] + dt*(Da1*x[1] - x[2] - Da2*x[2]*x[2] + u[2])
    return [x1 x2]
end

# Reactor 3
function f3(x::AbstractArray, u::AbstractArray)
    x1 = x[1] + dt*(-x[1] - Da1*x[1] + u[1])
    x2 = x[2] + dt*(Da1*x[1] - x[2] - Da2*x[2]*x[2] + u[2])
    return [x1 x2]
end

# general computation parameters
# d = 10

# # function compute_distributed(d::Int)
# Xc = [@interval(-5,5), @interval(-5,5)]
# d = 64
# dsc = DiscreteSpace(Xc,[d,d])

# Rkall = [[i,j] for i in 1:d for j in 1:d]

# 2d reactor state data
X1 = [@interval(0.0,1.0), @interval(0.0,1.0)]
ds1 = DiscreteSpace(X1,[d,d])

# 4d reactor state data
X2 = [@interval(0.0,1.0), @interval(0.0,1.0), @interval(0.0,1.0), @interval(0.0,1.0)]
ds2 = DiscreteSpace(X2,[d,d,d,d])

# 6d reactor state data
X3 = [@interval(0.0,1.0), @interval(0.0,1.0), @interval(0.0,1.0), @interval(0.0,1.0), @interval(0.0,1.0), @interval(0.0,1.0)]
ds3 = DiscreteSpace(X3,[d,d,d,d,d,d])
# computations

# reactor 1
u1 = [@interval(0.0,1.0)]
R1k = [[i,j] for i in 1:d for j in 1:d]

G1 = create_graph(R1k, ds1, u1, f1)
nlc1 = select_nonleaving_cells(G1)
lc1 = setdiff(collect(1:length(R1k)))
R1k = R1k[nlc1]
R1 = create_boxes(R1k,ds1)
writedlm("R1_cstr.txt",unique!(vcat([collect(vertices(b)) for b in R1]...)))

# reactor 2
u21 = unique!([i for (i,j) in R1k])
u21 = hull(ds1.dX[1][u21])

u22 = unique!([j for (i,j) in R1k])
u22 = hull(ds1.dX[2][u22])

u2 = [u21,u22]

R2k = [[i,j] for i in 1:d for j in 1:d]

G2 = create_graph(R2k, ds1, u2, f2)
nlc2 = select_nonleaving_cells(G2)
lc2 = setdiff(collect(1:length(R2k)))
R2k = R2k[nlc2]
R2 = create_boxes(R2k,ds1)
writedlm("R2_cstr.txt",unique!(vcat([collect(vertices(b)) for b in R2]...)))

# reactor 3
u31 = unique!([i for (i,j) in R2k])
u31 = hull(ds1.dX[1][u31])

u32 = unique!([j for (i,j) in R2k])
u32 = hull(ds1.dX[2][u32])

u3 = [u31,u32]

R3k = [[i,j] for i in 1:d for j in 1:d]

G3 = create_graph(R3k, ds1, u3, f3)
nlc3 = select_nonleaving_cells(G3)
lc3 = setdiff(collect(1:length(R3k)))
R3k = R3k[nlc3]
R3 = create_boxes(R3k,ds1)
writedlm("R3_cstr.txt",unique!(vcat([collect(vertices(b)) for b in R3]...)))

# subsystem 1 - R1+R2
function f12(x::AbstractArray, u::AbstractArray)
    xx1 = x[1:2]
    xx2 = x[3:4]

    x1 = f1(xx1,u)
    x2 = f2(xx2,xx1)

    return [x1[1] x1[2] x2[1] x2[2]]
end

u12 = [@interval(0.0,1.0)]
R12k = [[i,j,k,l] for (i,j) in R1k for (k,l) in R2k]

G12 = create_graph(R12k, ds2, u12, f12)
nlc12 = select_nonleaving_cells(G12)
lc12 = setdiff(collect(1:length(R12k)))
R12k = R12k[nlc12]
R12 = create_boxes(R12k,ds2)
writedlm("R12_cstr.txt",unique!(vcat([collect(vertices(b)) for b in R12]...)))

# subsystem 2 - R2+R3
function f23(x::AbstractArray, xi::AbstractArray, u::AbstractArray)
    xx1 = x[1:2]
    xx2 = x[3:4]

    x1 = f2(xx1,xi)
    x2 = f3(xx2,xx1)

    return [x1[1] x1[2] x2[1] x2[2]]
end

u12 = [@interval(0.0,1.0)]
R23k = unique!([[k,l,m,n] for (i,j,k,l) in R12k for (m,n) in R3k])

G23 = create_graph_4d(R23k, ds2, R12k, ds2, u12, f23)
nlc23 = select_nonleaving_cells(G23)
lc23 = setdiff(collect(1:length(R23k)))
R23k = R23k[nlc23]
R23 = create_boxes(R23k,ds2)
writedlm("R23_cstr.txt",unique!(vcat([collect(vertices(b)) for b in R23]...)))

# extract R3k
R13k = [[i,j,m,n] for (i,j,k1,l1) in R12k for (k2,l2,m,n) in R23k if k1==k2 && l1==l2]
R13 = create_boxes(R13k,ds2)
writedlm("R13_cstr.txt",unique!(vcat([collect(vertices(b)) for b in R13]...)))

# 3 reactors in series
# Rkall = [[i,j,k,l,m,n] for i in 1:d for j in 1:d for k in 1:d for l in 1:d for m in 1:d for n in 1:d]
# R1tt = select_test_cells(nlc12,lc12,G12,d,Rkall) 
R2tt = select_test_cells(nlc23,lc23,G23,d,Rkall)

# reconstruct sets
Rr = unique!([[i,j,k1,l1,m,n] for (i,j,k1,l1) in R12k for (k2,l2,m,n) in R23k if k1==k2 && l1==l2])

Rt = unique!([[i,j,k,l,m,n] for (i,j,k,l,m,n) in Rr if [k,l,m,n] in R2tt])
# validate sets
function f123(x::AbstractArray, u::AbstractArray)
    xx1 = x[1:2]
    xx2 = x[3:4]
    xx3 = x[5:6]

    x1 = f1(xx1,u)
    x2 = f2(xx2,xx1)
    x3 = f3(xx3,xx2)

    return [x1[1] x1[2] x2[1] x2[2] x3[1] x3[2]]
end

# Rt = Rr
Rkr = validate_set_6d(Rr, Rt, ds3, u1, f123)

# save data
bb = create_boxes(Rkr,ds3)
writedlm("R_cstr.txt",unique!(vcat([collect(vertices(b)) for b in bb]...)))
# Rt = unique!([[i,j,k] for (i,j,k) in Rr if [i,j] in R1tt || [j,k] in R2tt])

# end
println("Done!!!")

end

println("starting up")
main(4)

println("main computation begin")
@time main(16)

println("complete")

# end
# arr = ConvexHullArray(bxs)
# ak1 = Projection(arr,[1,2])
# ak2 = Projection(arr,[2,3])

# for plotting
# R1c = 